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Abstract. A new numerical method for solving a scalar ordinary differential 
equation with a given initial condition is introduced. The method is using 
a numerical integration procedure for an equivalent integral equation and is 
called in this paper an integrating method. Bound to specific constraints, the 
method returns an approximate solution assuredly within a given tolerance 
provided by a user. This makes it different from a large variety of single- and 
multi-step methods for solving initial value problems that provide results up 
to some undefined error in the form 0(h*), where h is a step size and k is 
concerned with the method's accuracy. Advantages and disadvantages of the 
method are presented. Some improvements in order to avoid the latter are 
also made. Numerical experiments support these theoretical results. 



1. An IDEA OF THE METHOD 

Consider a following differential equation supplied with an initial condition: 
f^=/(2/)-5(a;), O^x^b, 

\yio) = yo- 



The purpose is to find y{b), provided the solution can be extended from the ini- 
tial condition at a; = to a; = fe. Classical single or multi-step methods, e.g. 
Runge-Kutta or Adams schemes, would discretize the interval [0,6] into parti- 
tion {0,xi,X2t--,xn-i,xn} with xjv = b and compute iterative relations that 
would lead to some approximation of y{b). The answer would be in the form 
y(b) = yx^ + 0{hP), where y^^ denotes the approximate solution at point Xi of the 
partition, /i is a characteristic mesh size and p is the order of the method depending 
on the approximating scheme. The error 0{hP), in general, cannot be estimated. 
Some bounds may be found only in certain cases. These methods are mostly used 
as trustworthy, assumed p is large enough to drive the whole error 0{hP) to zero. 
Negative effects may happen if, for example, unstability takes place. 

Now, assume functions / and g from ([1]) satisfy the following conditions. 

(A) g{x)dx can be evaluated exactly. 

(B) Let t{z) = g{x)dx. Then t'(z) = g{z) > for any z > 0. 

(C) /(yo) > and f'{y) > for any y ^ yo- 
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Figure 1. Function ^(y) 



Let us call those the integrating conditions. From ([T]), using the separation of 
variables, it is easy to obtain an equivalent relation 

-— = / g{x)dx. 

yo fiy) Jo 

The term on the right-hand side is non other than t(&). By assumption (A) of 
the integrating conditions, this term is given precisely, i.e. without an error. From 
the numerical point of view, this and assumption (B) of the integrating conditions 
allow to replace r(6) with b without loss of generality to reduce problem ([T]) to 



y(o) = yo, 



(2) 

and 

(3) b. 



yo 



fiy) 



Denote 



p{y) 



fiy)' 



So ^ turns into 

rvib) 

(4) / p{y)dy = b. 



■' yo 

From now on, we will be focusing on integral equation (|1]). This is the core idea 
of the integrating method that will be presented below. The qualitative graph of 
p{y) is shown on picture [T] Such behavior is due to assumptions (C) and (D) of 
the integrating conditions. 

Remark 1. The condition that J^^°° p{y)dy = c is equivalent to the condition that 
the solution of ([2]) can only be extended to point x = c. Therefore, it is necessary 
to impose condition b < c, in order to deal with a solvable problem. 
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Figure 2. Lower rectangular and trapezoidal integration rules 



In order to solve @ for y{b), consider the lower rectangular and the trapezoidal 
methods with constant step h for the integral in (|H) , picture [H 

Denote the approximate integrals computed with the lower rectangular rule and 
the trapezoidal rule as J2i h n ^'^'^ '^t h n respectively, where N is the number of 
subintervals of length h. Formally, these sums are defined as 

N 

^ := ^h-p{yo + h-i), 

lJi,N i=l 
N ^ 

- II 2 ■ + P{ya + h-{i~ 1))). 

t,h,N i=l 

It should be clear that 

E = E +^>^^^' 

t,h,N l,h.N 

where Au^n = | • (p(yo) - p(yo + h ■ N)). 

Since function p{y) is concave up due to assumption (D) of the integrating con- 
ditions, the integration error on a single subintcrval [yo -\- h ■ i^y^) + h ■ [i + 1)] is 
simply an area of one of the regions presented on picture [3l depending on which of 
the two methods is currently used. 

More specifically, a region dounded by the blue lines and the black one corre- 
sponds to the lower rectangular rule and a region dounded by the green and the 
black lines corresponds to the trapezoidal rule. 

The following relation is obvious, but nevertheless is very important: 

(5) E < / p^y'^'^y < E +^h,N- 

l,h,N ''yo l,h,N 

It allows to formulate the idea of the integrating method as follows. For a given 
tolerance e find such natural numbers ni, n2, ni < n2, that 

• h ■ n2 — h ■ ni ^ e, 
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{yo + h- i, p{yo + h-i)) 




{yo + h- 1, p{yo + h ■ {i + 1))) {yo + h ■ {i + 1), p{yo + h ■ {i + 1))) 

Figure 3. Local integration errors for both rules 




yo + h-ni yo + h-n2 



Figure 4. The idea of the integrating method 

Once such pair is found, according to ([5]) we obtain 

/ P{y)dy <b < p{y)dy. 
J yo •'yo 

As J^^p{y)dy is an increasing function of z, it is clear that 

h ■ ni < y{b) < h ■ n2. 

Obviously, in either case 

\y{b) -h-ni\ < e, 
\y{b) -h-n2\ < e. 

Thus set the approximate solution of (U) as t/f, = h ■ ni or yf, — h ■ 112- For even 
better tolerance |, set the approximate solution ds yb — \ ■ h ■ (rii + 712), which is 
a mid-point between the previous two. 

In other words, we managed to find an approximate solution y of ((2) at point 
X = b within provided tolerance e, using the integrating method applied to 
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2. A STATEMENT OF TWO ALGORITHMS 

Once the idea of the method is clear, it is necessary to formulate an algorithm 
that realizes this idea consistently. Let superscript (j) correspond to the j-th iter- 
ation of the algorithm. 

Start with executing the lower rectangular integration of J^^^ p{y)dy with con- 
stant step h = h^^^ := e. Accumulate ^(i) until inequality ^(i) ^ ^ bis met. 
Denote such N as N = Then go one step back, i.e. to n^^ := — 1, and 
compute Et,Mi),„a)_i = Ei,Mi),n<^)-i +^fc(i),4^)-r ^^tM'Kn'^^ -1 ^ ^' then the 
approximate solution is found at the first iteration and yi, = • (2 ■ 712^' — 1). The 
algorithm terminates. However, it is possible that inequality E^ f^ii) „(i)_i ^ ^ does 

not hold. Thus it is necessary to try the second iteration. Set /i^^^ := = |. 
Redo the same lower rectangular integration and find Once again, it is im- 

portant to go back in order to check whcthcir the trapezoidal sum is less than or 

(2) 

equal to b. Only this time, make two backward steps instead of one, from to 
n'f^ := n-2'' ~ 2) as we halved h^^'' to obtain h^'^\ Depending on whether inequality 
Ef ?i(2) „<=^>_2 ^ ^ is satisfied or not, we either terminate the algorithm and set the 

approximate solution as yi, = • (2 • — 2) or continue to the next iteration, 

that is number 3. Note that for y^ we are returning a mid-point between the two 
neighboring ones for better accuracy |. It is not mandatory if we simply wish to 

reach tolerance e, and can choose either h'^'^ ■ n^^ or /i^'^ • n^'^ for y^. 
So the algorithm may be formalized as follows. 

Algorithm 1. At iteration j , execute: 

(1) hU) = i. 

(2) Find the smallest such integer n2, that E; ^(j) ^ b- Denote it as ri^^ . 

(3) If E/,M.),„0'-j+^/.0),„0'-j ^ ^' set y^ = .{2- nf - j) or y^ = 
h^^^ ■ {n2^ — j) or yb = /i^-^^ • n2^ . Either one returns an approximation of 
y(b) within tolerance e. Terminate the algorithm. 

(4) Else, continue with iteration j + 1. 

Theorem 1. The above algorithm has a finite number of iterations. In other words, 
it converges. 

Before we prove this theorem, it is necessary to prove the following lemma first. 
Lemma 1. 

for any iteration j . 

Proof. The lower rectangular integration rule applied to J p{y)dy creates subintcr- 
vals of length h^^^ = e each on y-axis, starting from y = yo. The integration at any 
of the following iterations j > 1 uses subintervals of length h^^^ = 4 each and thus 
contains nodes created by the first iteration in its set of nodes. It is also obvious 
that 

k-j 

i={k-l)-j+l 
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Figure 5. The lower rectangular integration for j ~ 1 and j = 3 

for any natural numbers k and j, since p{y) is a decreasing function. Picture [S] 
demonstrates this property in the case j = 3. From this, it follows that 

E ^ E ■ 

But 

E 

Since, by definition, is the smallest such natural number that 

E 



it is evident that 



j ■ 



From here, we obtain 

This implies the statement of the lemma. □ 

We are now ready to prove the theorem. 
Proof. Consider sequence aj — hU) n'-'^-j^ J > 1- It is evident that 
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Analyze these three terms separately. By definition of , inequality 

E <^ 

always holds. Also, according to Lemma 1 and the fact that p{y) is a decreasing 
function, the second term is bounded in a way presented below: 

^ • (p(2/o) - p(2/o + /i^'^ • n« - /iW)). 

The third term is non other than the lower rectangular sum with the step size 
h^^' between the nodes with numbers and rij — 1, taken with a negative sign. 
The following lower bound for it will be helpful. Picture [5] helps understand its 
meaning. 

^^^^ • p^yo + ^'''^ • *) ^ ^'(j'o + ^^^^ ■ (^2^^ - 1)) ■ ("2^^ - 1 - '4^^) • ^^^^ ■ 

As n^"'' = rij"'-' — j, we obtain 



Apply Lemma 1 and the fact that p{y) decreases: 

p(2/o + h^'^ ■ {v^i^ - 1)) ^ p(yo + /i^^^ • 

So we end up with 

a, < 6 + ^ • (p(yo) - p(2/o + h^'^ ■ - /i^))) - p(yo + /^'^^ • n^^ ^) ■ [j - 1) • h^^\ 

The algorithm terminates when inequality aj ^ 6 holds true. It will be then suffi- 
cient to require that 

^ • b(2/o) -p(2/o + h('^ ■ - /i(i))) ^ p{yo + /j(i) • nW) • (j - 1) • h^^l 
Isolate j to obtain 

.7. ■>! .1 ( p{y,)-p{yo + h(')-ni'^ ~hW) \ 

Thus, once j becomes large enough to satisfy ([7]), the trapezoidal sum satisfies 
inequality n'^^-j ^ ^ ^^"^ ^'^^ algorithm terminates immediately. □ 

Remark 2. For simplicity, condition ([7]) may be replaced by a stronger one: 
Of course, this may increase the cost. 
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Several observations must be pointed out. 

It is easy to see that the iterations in Algorithm 1 are independet in such a sense, 
that they can be carried out in any order since none of the iterations use the data 
obtained by the previous ones. They are presented in an order of increasing cost, 
since at every next iteration there are more nodes to evaluate p{y) at during the 
lower rectangular integration. 

Denote the smallest integer j satisfying ([7]) as js- The condition that j ^ js is 
sufficient but is not necessary for the convergence of the method. It may be that 
becomes less than or equal to b before j is large enough to satsify 
j ^ js- The smallest j of all those granting the convergence, denoted ja, is not 
known. This explains why Algorithm 1 starts with smaller j and does not jump 
right to js that assures the convergence with the highest cost of all. It seems 
this strategy may prove cost-saving if the right-hand side of (O is very large and 
therefore it makes sense to try fewer computations. However, as will be shown 
below. Algorithm 1 is only valuable from the theoretical point of view and must be 
avoided in real life computations due to availability of a better algorithm. 

Definition 1. is the largest such N , that satisfies 

Lemma 2. 

and 

Proof. By definition of rig"'-' , 
But 

E < E • 

Since 

E 

and 'Ylii h^i) TV is ^n. increasing function of iV, the first inequality immediately follows. 

To prove the second inequality, we use the idea similar to the one employed in 
the proof of Lemma 1. By topological properties of a trapezoidal sum applied to a 
function with a positive second derivative, it is clear that 

E < E 

for any natural number j. But 

E 
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(1) 

By definition of rig , we get 



Multiply by /i^^^ to obtain 



□ 



Theorem 2. For Algorithm 1 to terminate, it is necessary that j satisfies 
(9) J > 1 + J 



1 PiVo) - PiVo + h^'^ ■ - /i(i)) - 2 • pivo + h^^) ■ nW) 



Proof. Consider again the trapezoidal sum aj in ([B]) . The termination of Algorithm 
1 implies that aj ^ 6, or 

The right-hand side of this inequality is obviously bounded by ft-^^' ■p{yQ + h^^'' ■n'"/^), 
due to the definition of ■ Consequently, 

(10) ^;,o,.„o)_, - E ^^'^ • ^'(yo + ^^'^ • ^ f"^'^ ■ p(yo + h'^'^ ■ ^2^)- 

It is now necessary to find lower bounds for the two left-hand side terms. Interest- 
ingly, the proof of Theorem 1 would require to have zero on the right-hand side and 
the upper bounds of the left-hand side terms at their places. This would lead even- 
tually to ([7]), whereas we are working on a weaker condition that would apparently 
lead to dH). 

Clearly, due to Lemma 2, 



so, by decreasing behavior of ^(y). 



Next, by the same reason. 



h^'^ ■ P{yo + h^'^ ■ i) ^ p{yo + h^'^ ■ {ni''> + 1 - j)) ■ h^^^ ■ (j - 1) < 



for j > 1. Also, the right-hand side of ([T0| is bounded by • p{yo + h^^^ ■ rS^^). 
So from (flOl) there follows 



2j 



< 
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i 

Solving this inequality for j returns □ 



Remark 3. Condition ([9]) may be replaced by a stronger one: 
(11) j>l + 



1 p{yo)-3-p{yo + h(^^ -ni'^) 



2 2 ■ p{yo + hW ■ ni^^ - h(^)) ' 

Let the smallest integer j satisfying ([9]) be denoted jn ■ Having now the necessary 
condition of convergence provided by Theorem 2, it is easy to see that Algorithm 
1 is not efficient since it assumes than no information is provided on how fast the 
integrating method converges. Now that a lower bound for ja is known due to 
it will be wiser to focus on searching ja only among those j satisfying j„ ^ j ^ js- 
However, even knowing jn does not provide any better strategy of solving (|4]) than 
simply picking j := js immediately and doing the lower rectangular integration 
once ( not counting the one made with initial step h^^^ = e necessary for computing 
js via ^ ) with the maximum cost. This may be explained by the fact that the 
search for ja itself may result in a total cost higher than that of a single integration 
with the smallest step h^^-^K 

Suppose, for instance, the bisection method is used to find ja- At the first try, 
ji ~ i ■ (j„ + js)- Assuming that no information other than the values of j„ and 
js is provided, the probability that ji will grant convergence is about 50 %. The 
computational cost CbisecUon.i in this case is estimated roughly via the formula 

^ _ jn + js b b 

^bisection A ~ ^ ' ' ■ 

2 e e 

The last term comes from the very first integration with step h^^") — e. It is necessary 
in order to obtain j„ and js ■ At the same time, the computational cost Creai that 
corresponds to the case j = js is estimated via 

^ ■ b b 

^real Js ' ~\~ 

e e 

Obviously, '"^""''""■^ > 0.5, which means we do not even drop to the half of the 
maximum computational cost when we use j = ji . The probability of 50 % is not 
high enough to convince us to try ji. Even worse, already the second iteration 
j2 = 5 • (ii + js) with 75 % of success probability requires more computational time 
than the case with js since 

C bisection ,2 ~ ^ ' ~ ^~ Cbisection.l ^ Creal- 

So it is evident that the algorithm which uses as soon as js is provided, is 

the best choice from the point of view of computational efforts. It is formulated as 
follows. 

Algorithm 2. Perform the following two actions. 

(1) Execute the lower rectangular integration with step /i^^-' e until inequality 
hW n*^' ^ b is met. Check ifYlit ^{1) ^ b. If it is true, terminate 

the algorithm having set yb = h^^'' - n^2 ^ . Else, continue to step 2. 
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(2) Find the smallest such integer j that satisfies ^ or Compute /i^^-* = 4. 
Execute the lower rectangular integration the second and the last time, with 
step h^^\ setting j/f, — h^^^ ■ ri^'' . Terminate the algorithm. 

3. Advantages and disadvantages 

The main advantage of the integrating method (and apparently its right for 
existence) is that it returns solution y{\)) of ^ within provided tolerance e. This 
is what makes it unique and distinguished from the widely used Runge-Kutta or 
Adams methods that return the answer with an error in the form 0{h^). 

The other advantage is a stability of the method. Since it uses integration 
process, the method is always stable, regardless of given e and integration step 
h. This is not always true for other solvers of ordinary differential equations. For 
example, forward Euler method may experience unstability thus becoming unusable 
for large time-steps. 

Another obvious advantage is a simplicity of the algorithm. It is easy to under- 
stand and implement in any mathematical software. 

One of the disadvantages of the integrating method is that it is not general 
and only works for scalar problems having form ([1} and satifying the integrating 
conditions (A)-(D) imposed on functions /(y) and g{x). It is also important to 
know the value of Jy^ °° p{y)dy ( or at least some approximation of it ) in case this 
integral converges. It may be challenging to use the integrating method without 
being sure whether the problem is solvable in general or solvable within a reasonable 
amount of time. The latter may be an issue if b is very close to c which is the right 
limit of the solution extension interval [0,c). A relatively high cost of the method 
comes from the integrating process with small step size and is considered another 
disadvantage. 

Remark 4. All these advantages and disadvantages remind of similar situation 
with numerical methods for solving nonlinear algebraic equations. The bisection 
method is known for its ability to return the solution of f{x) — on interval [a, b] 
within a given tolerance e as long as f{a) ■ f{b) < 0. The method, however, cannot 
be extended to a class of vector equations. It also does not converge as fast as the 
Newton's method. The latter can also be used for vector equations. But in return, 
it cannot assuredly present the answer within the given tolerance e. 

Another property of the integrating method that may seem as a disadvantage 
is that it works towards finding y{b) only, whereas all known step solvers return 
intermidiate points as well. They therefore present an approximation of the whole 
solution curve y{x) on the presented mesh Xk, < Xk ^ b. The fact that the inte- 
grating method manages to find y{b) within tolerance e does not necessarily imply 
that intermidiate points y{xk) are approximated within e too. In order to compute 
these intermidiate points, the integrating method may be applied separately to each 
point. This will result in a very high total cost. 

The last problem may be avoided with help of inequality ([7]). 

Theorem 3. Let {0, cci, X2, xat-i, xat} with xn — b be the mesh provided by a 
user and n^^^ is the smallest such integer that ^^(^j ^(i) ^ b with h^^^ — e. Let 
j be an integer that satisfies ([7]) or ([S]), and the lower rectangular integration is 
executed once with step h^'^^ = |. For every Xi set y^^. ■= h^-^^ ■ Hj"'^, where n^], is 
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the smallest integer such that ^u) ^ ^^^^"^ following holds true: 
\y{xk)-y^^,\ < e Vfc, 1 ^ fc iV. 

The above theorem assures that if j satisfies (O or ([8]), concerned only with 
the final node xn — b, the solution at the intermidiate mesh nodes is already to 
be evaluated within tolerance e during the lower rectangular integration aiming to 
find yb- This is a very relieving result. 

Proof. Pick any intermidiate point Xk, < Xk < b. When the integrating method 
is used specifically to find y^;^, then it will be sufficient to use the lower rectangu- 
lar integration with step size h^^''^ — j^, where integer jk is the smallest integer 
satisfying 

. ^ ^ if p{yo) - Pivo + e • "2^1- - e) \ 

Jk^^+n-[ : mT ' 

^ \ P(yo + e-"2,fc) / 

with 71^:^1 being the smallest such integer that satisfies J^i ^ ^ ^k- Since Xk < b, 
it is obvious that n':^l ^ nj^]^ — n^2 ^. So 

PiVa) - Pivo + e • '4^^ - e) ^ Pivo) ^ PiVo + e • "2!! ~ 
By the same reason, 

p{yo) - pivo + e • 4^' - e) pivo) - PiVo + e • - e) 



So we see that if j satisfies it automatically satisfies 

^. ^ ^ ^ 1 / p(yo)-p(?/o + £-"2^fc-£) \ 

^ V P(yo + e-4!fe) / 
This proves that j ^ jk- This assures the integration step /i*^-'^ = is small enough 
to guarantee that \y{xk) — yx^l < ^ fo^' arbitrary Xk, < Xk < b. □ 

Remark 5. Similar proof may he conducted in case of inequality ([8]). 

Remark 6. Similar result exists for the necessary condition It simply states 
that the larger integer k is, the stronger condition ([9]) for point Xk is. In other 
words, if ([9]) is satisfied for the final point h, it is automatically satisfied for the 
previous points Xk of the mesh. The same may he said in case of inequality (jlip. 
This result seem to present no value for the current work though. 

Now that Theorem 2 is proven, Algorithm 2 may be assuredly used to find the 
approximate solution of ([2]) on arbitrary mesh {0,xi, ...,xn-i,xn}, xn = b, with 
almost the same cost as for solely the final point b. 

Conjecture 1. The integrating method may he improved to cover a wider area of 
initial-value problems than just those restricted hy the integrating conditions (A)- 
(D). If it is possible to track all the inflection points of function from ([1]) 
precisely, it may he possible to use a lower rectangular integration on intervals 
where 

1 \ " „ 
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and a higher rectangular integration on intervals where 




< 0. 



Using relations analogous to ^ may help construct algorithms to find solutions 
within tolerance e even for this more general class of problems. This may he a topic 
of further research. 



4. Numerical experiments 

Numerical experiments were conducted to test both Algorithm 1 and 2 on two 
initial- value problems: 



(12) 



[2/(0) = 0, 

with exact solution y{x) — — 1 and mesh Xk — 0.05 ■ k, k ^ 20, and 



(13) 



dy 

y(0) = 0.5, 



with exact solution y — and mesh — 0.05 • fc, fc 5C 32. It is important to 
note that in the second case the solution may only be extended up to a; = 2 since 
the integral 



Jo.5 

converges and is equal to 2. In the first case. 



+ 00 



dy 

= +0O, 



y + 1 



so no problems are encountered. 

Tolerance was e = lO"**. Both algorithms were implemented in the mathemat- 
ical software Octave. Algorithm 1 was applied separately to each node Xk- The 
computational time was measured on a computer with Intel Core 17 processor with 
2.80 GHz frequency. The computation started with doubled tolerance h'-^^ — 2e 
and the approximate solution y^^ was found as a mid-point between yo + h^^^ ■ rt^^'-* 
and yo -t- /i'--'-' • in both cases. Tables 1 and 2 present results of Algorithm 1 
applied to problems (fT2|) and (|13p . The first case required 4.3 seconds, and the 
second one required 75.6 seconds. 



14 ALEXANDER V. LOZOVSKIY 







Iz/fxt) — \ • 10"* 

1 y y-^K. } \ 


The actual number of iterations 


jn 


J s 


0.05 


0.0513 


0.289 


1 


1 


2 


0.1 


0.1051 


0.709 


1 


1 


2 


0.15 


0.1619 


0.658 


1 


1 


2 


0.2 


0.2215 


0.972 


1 


1 


2 


0.25 


0.2841 


0.746 


1 


1 


2 


0.3 


0.3499 


0.412 


1 


1 


2 


0.35 


0.4191 


0.326 


1 


1 


2 


0.4 


0.4919 


0.753 


1 


1 


2 


0.45 


0.5683 


0.122 


1 


1 


2 


0.5 


0.6487 


0.213 


1 


1 


2 


0.55 


0.7333 


0.4698 


1 


1 


2 


0.6 


0.8221 


0.188 


2 


1 


2 


0.65 


0.9155 


0.408 


2 


1 


2 


0.7 


1.0138 


0.473 


2 


1 


2 


0.75 


1.1171 


0.9998 


1 


1 


2 


0.8 


1.2256 


0.591 


2 


1 


2 


0.85 


1.3397 


0.531 


1 


1 


2 


0.9 


1.4597 


0.969 


1 


1 


2 


0.95 


1.5857 


0.097 


2 


1 


2 


1.0 


1.7183 


0.182 


2 


1 


2 



Table 1: results of Algorithm 1 applied to each Xk for solving ([T2l) 



Results of Algorithm 2 applied to p2|) and (fT3)) are presented in Tables 3 and 
4 respectively. The initial step /i'-^^ = 2e was divided by 2 in case (fT2|) and by 14 
in case This agrees with the cell value of the last line and the last column of 

Tables 1 and 2 respectively. Algorithm 2 significantly speeds up computations since 
it only took 0.96 seconds to obtain the results of Table 3 and 6.04 seconds for Table 
4. This may be explained by the fact that, although inequality ([T]) is not a criterion 
of the algorithm termination, it returns a value close to the actual minimum number 
J by which h^^^ must be divided in order to grant required tolerance at every fixed 
node Xk ■ Both Tables 1 and 2 show that the difference between and js is no larger 
than 1 for all nodes. If the fact that j„ and js are close to each other remains true in 
many other cases, it is another plus towards Algorithm 2. The author recommends 
using it in real life applications where a need to reach the desired accuracy with 
100 % guarantee is higher than a need for computational speed and efficiency. 
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V-r, 


\v(xk) ~ Vt, I ■ 10^ 


The actual number of iterations 


7n 


?\ 
J ^ 


0.05 


0.5129 


0.795 


1 


1 


2 


0.1 


0.5263 


0.158 


1 


1 


2 


0.15 


0.5405 


0.405 


1 


1 


2 


0.2 


0.5555 


0.556 


1 


1 


2 


0.25 


0.5715 


0.714 


1 


1 


2 


0.3 


0.5883 


0.647 


1 


1 


2 


0.35 


0.6061 


0.394 


1 


1 


2 


0.4 


0.6250 


0.000 


2 


1 


2 


0.45 


0.6451 


0.613 


2 


1 


2 


0.5 


0.6667 


0.333 


1 


1 


2 


0.55 


0.6897 


0.448 


1 


1 


2 


0.6 


0.7143 


0.143 


1 


1 


2 


0.65 


0.7408 


0.593 


2 


1 


2 


0.7 


0.7693 


0.692 


1 


1 


2 


0.75 


0.8000 


0.000 


2 


1 


2 


0.8 


0.8334 


0.667 


2 


1 


2 


0.85 


0.8696 


0.348 


2 


2 


3 


0.9 


0.9091 


0.091 


3 


2 


3 


0.95 


0.9524 


0.524 


3 


2 


3 


1.0 


1.0000 


0.333 


3 


2 


3 


1.05 


1.0527 


0.684 


3 


2 


3 


1.1 


1.1112 


0.556 


3 


2 


3 


1.15 


1.1766 


0.961 


3 


3 


4 


1.2 


1.2501 


0.500 


4 


3 


4 


1.25 


1.3334 


0.667 


4 


4 


5 


1.3 


1.4286 


0.486 


5 


4 


5 


1.35 


1.5385 


0.785 


5 


5 


6 


1.4 


1.6667 


0.619 


7 


6 


7 


1.45 


1.8183 


0.896 


7 


7 


8 


1.5 


2.0001 


0.778 


9 


8 


9 


1.55 


2.2223 


0.978 


10 


10 


11 


1.6 


2.5001 


0.857 


14 


13 


14 



Table 2: results of Algorithm 1 applied to each Xk for solving ([T^ 
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Vxk 


- • 10^ 


0.05 


0.05120 


0.711 


0.1 


0.1051 


0.709 


0.15 


0.1618 


0.342 


0.2 


0.2214 


0.028 


0.25 


0.2840 


0.254 


0.3 


0.3498 


0.588 


0.35 


0.4190 


0.675 


0.4 


0.4918 


0.247 


0.45 


0.5683 


0.122 


0.5 


0.6487 


0.213 


0.55 


0.7332 


0.530 


0.6 


0.8221 


0.188 


0.65 


0.9155 


0.408 


0.7 


1.0138 


0.473 


0.75 


1.1170 


0.0002 


0.8 


1.2256 


0.591 


0.85 


1.3397 


0.531 


0.9 


1.4596 


0.031 


0.95 


1.5857 


0.097 


1.0 


1.7183 


0.182 



Table 3: results of Algorithm 2 applied to 
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\y(xi.) -Vt\- lO'^ 


0.05 


0.5127 


0.919 


0.1 


0.5262 


0.872 


0.15 


0.5404 


0.977 


0.2 


0.5555 


0.841 


0.25 


0.5713 


0.8571 


0.3 


0.5881 


0.924 


0.35 


0.6060 


0.892 


0.4 


0.6249 


0.857 


0.45 


0.6451 


0.899 


0.5 


0.6666 


0.810 


0.55 


0.6896 


0.837 


0.6 


0.7142 


0.857 


0.65 


0.7407 


0.836 


0.7 


0.7691 


0.879 


0.75 


0.7999 


0.857 


0.8 


0.8333 


0.762 


0.85 


0.8695 


0.795 


0.9 


0.9090 


0.766 


0.95 


0.9523 


0.810 


1.0 


0.9999 


0.714 


1.05 


1.0526 


0.744 


1.1 


1.1110 


0.683 


1.15 


1.1764 


0.563 


1.2 


1.2499 


0.571 


1.25 


1.3333 


0.476 


1.3 


1.4285 


0.429 


1.35 


1.5384 


0.330 


1.4 


1.6666 


0.238 


1.45 


1.8182 


0.104 


1.5 


2.0000 


0.143 


1.55 


2.2223 


0.349 


1.6 


2.5001 


0.857 



Table 4: results of Algorithm 2 applied to (fT3|) 
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